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PACS numbers: 02.70.-c, 03.65.-w, 03.75.Lm, 34.20.Cf 



* Electronic address: |rdecay@googlemail.com 
^Electronic address: | J . Br and@ massey. ac . nz 



1 



I. INTRODUCTION 



The decay of a particle by tunneling through a potential barrier into a continuum is 
a fundamental and unique phenomenon in quantum mechanics. The tunneling of multi- 
particle systems is just as important and presents one of the places where the understanding 
of macroscopic quantum phenomena can start [l| . The tunneling and decay of Bose-Einstein 
condensates (BECs) are attractive subjects of study {2], since the BEC is a unique state of 
matter where quantum mechanical features are manifested at the macroscopic level. After 
BECs were first realized experimentally in dilute atomic gases jsl, a huge amount of related 
research followed. Ultra-cold atoms are usually trapped in a finite potential well and the 
decay by tunneling into a continuum is an existing and potentially desirable possibility. In 
this context it was realized that understanding the decay dynamics by tunneling is a very 
important task 

In most cases BECs have thousands to millions of particles and the dynamics is ade- 
quately described by the Gross-Pitaevskii (GP) equation of mean-field theory [8|, a 
nonlinear Schrodinger equation. The GP equation governs the time evolution of phase and 
particle number density of an essentially fully Bose-condensed system. With many works 
on the mean- field description of BEC tunneling jo-lS], it is remarkable that there is still 



a discussion, both about the technical implementation 12| and the correct formulation of 

n 

mean- field theory related to the decay problem [10| . It is thus desirable to obtain a detailed 
understanding of the microscopic physics of multi-particle decay. 

The cases of stronger interactions or fewer particle numbers are also important, where 
the GP equation is less accurate. In the few boson regime, correlated decay of particles was 
observed and studied both experimentally and theoretically fl^ . [l5| . The particle correlation 
in the decayed wave is important in relation to the atom laser [l6|. For strongly- interacting 
bosons in a one-dimensional trap Bose-condensation is not relevant but the gas acquires 
properties related to fermionic systems In the Tonks-Girardeau limit of infinite inter- 
actions the few boson decay problem was treated analytically [l8| , and numerical simulation 
have addressed the crossover for finite interactions from a harmonic trap with up to four 
bosons 19|. The analytic treatment of few boson decay with finite interaction strength 
remains a difficult task. 

In this paper, we approach this problem by both numerically and analytically. We study 



2 



the simplest case of two repulsively interacting bosons in a potential trap in one dimen- 
sion. The time evolution of the decay is obtained from first principles by solving the time- 
dependent Schrodinger equation numerically. Then, it is compared to approximate analytic 
methods, starting from the exact solutions of local spatial regions. The decay phenomena 
are investigated for a wide range of interaction strength, from zero to very strong repulsion. 
The analytic model predicts exponential decay mode of the interacting system, which is in 
very good agreement with our numerical simulation. Also, the decay of the total particle 
number is well explained with a simple theoretical model. 



II. THE MODEL HAMILTONIAN 



We choose a model Hamiltonian for the two interacting boson decay. Considering the 
kinetic energy, external potential Vg^ for trapping and interaction U between bosons, the 
total Hamiltonian with two identical bosons is written as 

The external potential is 

{oo for a; < 

(2) 
V6{x -L) for 5 > 

and acts as potential trap by a delta barrier at position L. This choice of external potential 
has some advantages in that the delta barrier has zero width so the consideration of decay 
process inside the barrier is not needed. Also, the analytical treatment of the decay process is 
simplified. In single particle case we found that the exponential decay mode dominates and 
non-exponential features are strongly suppressed compared to a finite-width barrier case. 
Computationally, the narrow width of the delta function makes the Hamiltonian matrix 
more sparse, which makes the problem tractable. 



Considering only s-wave scattering 
particle 2 is given as 



3, 



the interaction potential between particle 1 and 



U{xi,X2) = g5{xi- X2) (3) 
where ^ is a coupling constant and xi and X2 are the positions of each boson, respectively. 
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To simplify the analysis and compare the result with external parameters, we introduce 
dimensionless units. The new length unit x is defined a.s x = x/L. The Hamiltonian is 
rewritten as 

+ - X2). for a; > (4) 

Dividing both sides by f?j{mL'^\ we get the rescaled, dimensionless Hamiltonian H = 

1 92 1 ^2 

H = --77^ - + V6{xi - 1) + V6{x2 - 1) + gS{xi - X2) (5) 



2 dxi 2 dxl 



Here 



mL ~ mL , . 

The Schrodinger equation with this Hamiltonian is given by 

id^^ = Hij, (7) 

where t = ht/{mL'^) with t is unsealed time. 

In Xi — X2 space, the Hamiltonian looks like figure [1] The dotted lines represent delta 
potentials from trap and interaction between particles. From now on we denote the region 
where both particles inside the trap as region (1), where one particle in and one particle out 
of the trap as region (2), and both particles out of trap as region (3). 



III. NUMERICAL SIMULATION OF TWO BOSON DECAY 

Now, we set up the decay of two interacting identical boson in this Hamiltonian. We 
choose the initial condition that both particles inside the delta trap as the two boson ground 

n 

state of y — 7- 00 case. Specifically, this initial state 'j/'inila^i, 2^2) is given by [20| 

^ini(a;i,X2) = Ni^i{{Ai{ku, k2i)e''''^''' - ^(fcH, k2i)e-''''^''') sm{k2iX2) 
HMku, fc2.)e*'""^ - A,iku, A;2.)e-*'=-^'^) sin(A;i,X2)) 
for < 0:2 < < 1 (8) 
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with 



A^iki, k2) 
kii and k2i satisfy the equation 



= {iki + ik2 + g){iki - ik2 + g) 
-- {iki - ik2 - g){iki + ik2 - g) 
-- {iki + ik2 + g){iki - ik2 - g) 
{iki + 'ik2 - g){iki - ik2 + g). 



kit = + arctan(- — — — ) + arctan( ^ 



'k 



k2i = — arctan( 



li — rv2i 



kli — k2i' 



+ arctan( 



kli + k2i 



kli + k2i ' 



(9) 
(10) 

(11) 
(12) 



(13) 
(14) 
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X, 



FIG. 1: Hamiltonian in Xi — X2 space. 
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The initial wave function in < xi < X2 < 1 region is obtained from the boson symmetry 
condition '?/'ini(a;i, 2:2) = V'ini(a;2, a^i)- In other regions the initial wave function is zero. The 
normalization constant A^ini is chosen to satisfy J dxidx2\ipini\'^ = 1- ^li and versus 
interaction strength g is shown in figure [2J For zero interaction both kn and are same 
as 7r, the single particle ground state wavevector. For nonzero g they rapidly deviate form 
n as g increases, and ku approaches to 271 and k2i approaches to vr (Figure [2]) . 




20 40 60 80 100 120 



g 

FIG. 2: kii and k2i versus g plot. The dashed lines are vr and 27r, the wavevectors of 
V = 00 lowest and next lowest states. As the interaction strength g increases, ku 
approaches to 2tt and k2i approaches to vr. 



To analyze the decay of interacting bosons, we solve the Schrodinger equation directly. 
The Schrodinger equation and its formal solution are 

idti^ = Hil) (15) 
V^(t) =exp(-iift)^(0). (16) 



We use Crank-Nicolson method to solve this equation numerically |2lL |22|. 

For the numerical representation of Hamiltonian, we choose the triangular region < 
X2 < xi < Xmax in X space (0 < xi < X2 < Xmax region is determined due to the bosonic 
symmetry), with X^ax is large enough that in our observing time very little decay products 
reach near X^ax- This region is discretized by dividing X^ax by N^, and all points in the 
triangular region are arranged in one column vector. The Hamiltonian matrix obtained by 
discretization of x space and using a finite-difference formula for the second derivatives can 
be quite large, but it is a sparse matrix as most elements are zero. 
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For small dt, 



exp(iM/2) ^(t + dt) = exp{-iHdt/2) tpi^t) 

ilj{t + dt) = (1 + iHdt/2)-\l - iHdt/2) ^(t) + 0{dt^). 



(17) 
(18) 



This method is second order in dt and unitary (i.e. probability is conserved). This is an 
implicit method, since it contains the inverse operator. The matrix inversion is efficiently 
implemented by solving linear equations. The time evolution of the wave function is obtained 
by iterating equation (|T8|) . 

For the simulations in next sections, the following parameters are used. Xmax = 45, 
Ax = (Xmax/^x) = 1/24, dt = 0.002, V = 5. The convergence of the numerical solutions 
is checked by changing spatial grid size and time step. We also check numerical simulation 
with known analytic solutions for special cases g = and g = oo. To see the effect from 
the reflection of waves at the boundary the results are examined by changing X^ax and by 
putting absorbing potentials near Xmax- In our parameter regime, those effects are very 
small and do not change the main results. 

IV. RESULTS AND ANALYSIS 

For the understanding of the decay of two interacting bosons, a good starting point is the 
parameter region where we know the exact analytic solutions. In our case, we know exact 
eigenf unctions of Hamiltonian for two extreme cases, g = and g = oo. In those cases, the 
two particle eigenfunctions are obtained by the combination of one particle eigenfunctions, 
which are known in analytic form. For arbitrary g > 0, the results lie between these two 
extremes, and the exact analytic forms are not known. 

A. Vanishing and infinite interaction limits 

When g = there is no interaction between two particles. They act independently, with 
only a symmetric wavefunction condition. The eigenfunction is written as 



ij{ki,k2,Xi,X2) 



1 



{(f){ki, Xi)(p{k2, X2) + 0(/c2, xi)(f){ki, X2)), 



(19) 



V2 
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where the total eigenenergy is £^ = {kj + A;|)/2 and (f){k, x) is the one particle eigenfunction 
with eigen wavevector k. In our model the explicit form of is given by 

{cAk) sm(kx) for < ,T < 1 

(20) 
C2(/c)e"^^ + C3{k)e for 1 < x 

where 



ci{k) = d- ^ , (21) 

Wl + ^ sin /c cos k + ^5- sin^ k 

C2(^) = ^(-(^ + ^) + ^e-2^')ci(A:), (22) 

cs{k) = l(^{t-j) + je'^'ym- (23) 

The one particle decay rate can be calculated by outgoing boundary condition, setting the 
coefficient of outgoing wave cs{k) — and solving for k (this is also the pole of scattering 
matrix). The equation Cs{k) — has complex solutions, each of them corresponds to different 
decay modes. We denote the complex solution of C3{k) = as kzo, kzi, ■■■ with kzo the lowest 
decay mode and kzi next lowest decay mode, etc. For the V — 00 ground state initial 
condition 

tpini{x) = v^sin(7rx), (24) 

the dominant decay mode is kzo- Since the decay mode wavefunction is also a complex 
eigenfunction, its time dependence is given by e~*^*, where E — A;^q/2. The time evolution 
of one particle probability inside the potential trap Piin(^) follows the exponential decay 

Piin{t) ~ |e-*^f = e-^^°^ (25) 
7^0 = -'^kzorkzoi, (26) 

where kzor and kzoi are the real and imaginary parts of kzo, respectively. 

The decay of two interacting bosons in the special cases of = and = 00 is obtained 
from the single particle decay patterns, respectively. 

For g — 0, the two particle wave function is the product of one-particle wave functions, 
and their decay is just the product of the individual decay. With the condition that the 
initial wave function was the ground state of ^ = 00: 

■0ini(a^i,a;2) = 2sin(7rxi) sin(7rx2). (27) 



If we write the probability of both particle inside the trap as Pi, probability of one particle 
in and one out as P2 and both particles out as P3, their dominant time evolutions are 

Pi(t) ^ e-2^^»*,, (28) 
PaW ~2e"^^«*(l-e-^^°*), (29) 
Psit) ^ (1 - e-^^"*)2. (30) 

Another case we know the exact eigenfunction of Hamiltonian is = 00 case. In this 
case, the two particle eigenfunction is written as 

Ipi^ki, k2, Xi,X2) = -^{(f){kiXi)(j){k2X2) - (t>{k2Xi)(t>{kiX2)), 

v2 

for xi > X2, (31) 

and ^/'(fci, ^2, xi, X2) = tp{ki,k2,X2,xi) for the Xi < X2 region. Like in the case of fermions 
the probability density is zero along xi = X2 line. The Vq = 00 ground state initial condition 
is given by 

'?/'ini(a;i, X2) = -\/2(sin(7rxi) sin(27rx2) — sin(27rxi) sin(7rx2)) 

for xi > X2 (32) 

and '?/'ini(a;2, a;i) = '?Aini(a;i, X2) for xi < X2- The decay oi g = 00 two bosons involves two 
decay mode, with lowest wavevector kzo and next lowest one kzi- The time evolutions of Pi, 
P2 and P3 are 

Pi(t) ^ e-(''^o+^^^)* (33) 
P2{t) ~ e-^^«*(l - e-"'^'') + e-^^^*(l - e"^^"*) (34) 
Psit) ^{1- e-^^''*)(l - e-^^i*) (35) 

with 

7zj = ~'2kzjrkzji (^^) 

and kzjr and kzji are real and imaginary parts of kzj, respectively. 
B. Arbitrary g > case 



For the general case of < < 00 exact analytic eigenfunctions are not known. We 
use the numerical method of section [Till to obtain the decay of probabilities. To conduct 
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the simulation, first the initial condition was chosen as the ground state of trap potential 
strength V = oo limit. 

Quite interestingly, the numerical results in this section show that a rather simple model 
can be used to explain interacting boson decay. For the decay of interacting particles, it is 
expected that the number density of particles shows non-exponential decay. When there are 
more particles inside the trap it decays faster, and with less particles the decay is slower. 
But if we examine the probability Pi of both particles inside and the probability P2 of one 
particle inside and another out separately, they show quite distinctive features. 

If we plot the logarithm In Pi vs time for various interaction strength g, the graphs 
show straight lines, meaning the decay is exponential. Furthermore, the decay rate can 
be obtained by theoretical estimation. Like the decay rate calculation of one-particle case, 
we can apply the outgoing boundary condition for the wavefunction in region (1). Since 
the probability of both particles escaping simultaneously is very small due to the repulsive 
interaction, it is ignored. Then the outgoing boundary condition from region (1) to region 
(2) can be written as follows. 

First, the wavefunction in region (1) satisfying the Bethe ansatz and the boundary 
conditions at Xi = and Xi = X2, can be written as [the form of the coefficients without 
normalization is given in equations (Q to ( [T2|) ] 

V'{i)(xi,X2) = iAiih,k2)e''^''' - A2{h,h)e-'''''')sm{hx2) 
+ iAsih, k2)e"''''' - A^ih, k2)e-'^'''')sm{kiX2), 

for < X2 < xi < 1 (37) 

and the outgoing eigenfunction in region (2), 1^(2)7 can be written as 

^(2)(xi,X2) = Pie'^i"^ sin(A;2a;2) + B2e"'''^' sin(fciX2), 

for 1 < xi, < X2 < 1 (38) 

with the boundary condition 

^(i)(l,X2) = V'(2)(l,a;2), (39) 
5x1^(2) (1,2:2) - <9^i^(i)(l,a;2) = 2V"^(i)(l,X2). (40) 

Conditions (l39ll and (HOil yield four equations with four unknowns Bi, B2, ki and /c2- Solving 
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for ki and k2 we get two equations 

A,ik,, k2)e'''' - A^ik,, A;2)e-*^^ = --yA^ik,, k2)e-''' (41) 

As{ki, k2)e''' - A,{k,, k2)e-'''' = -'-yA.ik,, k2)e-'^\ (42) 

and two complex wavevectors kig and k2g for their solutions. When we write real and imag- 
inary parts of complex eigenvectors as kig = kig^ + ikigi and k2g = k2gr + ik2gi, both of their 
imaginary parts are negative. Considering that the time evolution of an energy eigenfunction 
follows e~*'^* like the one particle decay mode, it can be expected that exp(— ^(A;^^ + fc|g)t/2) 
dominates in time evolution. When we compare the probability of both particle inside 
Pi{t) with \exp{-i{kjg + klg)t/2)\^, indeed we see that this is what happens. Both are in 
very good agreements as shown in figure [3l -Pi(t) decays exponentially with the decay rate 
predicted by outgoing boundary conditions. 

Pi(t) ^ I expi-tikl + kl)t/2)\' = e-^^\ (43) 
7g = -2kigrkigi " 2 A;23r j • (44) 

Figure H] shows 7^ change for various g. 7^ changes a lot for small g, and approaches to 
the decay rate of = 00 slowly. Figure [3] shows the comparison between —jgt line from 
theoretical prediction and log Pi from numerical simulation. They match very well well for 
all g > 0, thus showing Pi decays exponentially even with interaction between bosons. 

t 

-0.5 

In P, -1.0 
-1.5 
-2.0 
-2.5 
-3.0 
-3.5 

FIG. 3: InPi(t) plots for various g (different colors). Dots are from numerical simulation 
and lines are from theoretical prediction of decay rate by outgoing boundary conditions of 
( 14T1) and f H2|) . The numerical simulation and theoretical prediction show very good 

agreement. 

Next we consider the time evolution of P2, one particle in and one particle out of trap 
probability. It is more complicated than that of Pi, since it contains probability inflow from 
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region (1) and outflow into region (3). Like Pi case, we already know the dominant parts of 
P2{t) for special cases, g = and g = oo. 
For g = 0, the decay of P2{t) has the form 



P2,g=o{t) ^ 2e-^^«*(l - e 



-7zo^^ 



and for g = oo, 



2,9=00 



it) 



(45) 



(46) 



where 7^0, 7zi are the lowest and next lowest decay rate of one particle in the potential trap. 
For the g = case, both bosons decay from the same mode independently. For the g = 00 
case, two bosons decay from the separate decay modes without interfering since they are 
almost orthogonal. For general < g < 00, the time evolution of P2{t) will be between (H5|) 
and fH6l) and as g is increased P2{t) will change from (145|) to (1461) . We try different models 
for two regimes where g is not large (weak or moderate repulsion) and where g is very large 
(strong repulsion), and investigate regions of validity for each model. 

For the weak or moderate repulsive interaction, we try a simple model for P2 decay. If 
we assume that the probability of both particle escaping simultaneously is very small, which 
is satisfied when the decay rate is small and interparticle interaction is repulsive, then the 
inflow from region (1) has very simple form since the dominant part of Pi satisfies (H3|) and 
almost all escaping probability from region (1) goes to region (2). We can write P2 as 

dP2 r. 



dt 



out 



(47) 



where Pin is the probability inflow from region (1) to region (2) and Pout is the probability 
outflow from region (2) to region (3). 



Rate 




FIG. 4: Pi decay rate 7^ vs g plots. Solid line represents 7^, lower and upper dashed lines 
represent Pi decay rates of g = and g = 00 cases, respectively. 
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Fin is simply 7ge~'''s*, which is -Pi(t) outflow from region (1). For the form of outflow -Pout, 
we try exponential decay model. In that case, Fout is set as — 723-P2 where 723 is the decay 
constant from region (2) to region (3). With this assumption, the solution of ( H7|l has the 
form 



The decay constant 723 is yet undetermined, so (HSj) becomes one parameter fltting model. 
This exponential decay model of Font implies that remaining particle in the trap will decay 
exponentially like one particle decay after other one escapes, with only one decay mode. 

Compared with numerical simulation, model fHHl) shows very good agreements. Further- 
more, it shows that even for larger g the fltted parameter 723 is very close to the lowest single 
particle decay rate 720- Figure shows the comparison between numerical simulation and 
(HSj) with 723 substituted by 720 (dashed red line) and fHSjl with 723 obtained from fltting 
(blue line), for (7 = 0, 1, 10. All shows very good agreements and blue lines are not shown 
well due to overlapping with red. The agreements for even (7 = 10 is quite surprising, since 
for (7 = 10 the initial wavevectors inside trap are far from lowest decay modes as shown 
in flgure [2l The initial two wavevector ku = 5.347 and = 2.720, the escaping complex 
eigenvectors kig = 4.996 — 0.1445i and /c2g = 2.507 — 0.04881i (up to 4 signiflcant digits) 
for g = 10. kii and kig are closer to second decay modes, but still P2 decay to region (3) 
is dominated by single particle lowest decay rate. Figure [6] shows 723 compared to 7^0 and 
their relative differences for various g. It shows that the relative difference between 723 and 
7x0 are less than 2% for < g < 17, and the difference increases and approaches to 10% for 
larger g. 

In the strongly repulsive interaction region where the difference between 723 and 7^0 
increases, the deviation of model HH] from numerical simulation also increases. In this region 
we try different model which is close to (146|) . Physical meaning of (l46l) is that there are two 
decay modes and two decay modes decay independently without interfering each other. In 
our case, we have two complex eigenvector kig and k2g from (HTj) and (H2l) . Assuming that 
P2 decays from each complex wavevector and each mode do not interfere each other, we 
write decay model of P2(t) for large g as 



(48) 



Ig - 723 



(49) 
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where 

7lg = —'^kigrkigi, 72g = —2k2grk2gi (50) 

This model works better for larger g than fHSl) as figure [7] shows. The fH9|) model describes 
the peak of P2{t) well, and discrepancy with the numerical data becomes smaller for larger 
9- 

To see the agreements between numerical simulation and fitting model quantitatively, 
we consider the absolute mean of relative error rj between the probability calculated by 
numerical simulation P2num and the probability calculated by model P2modeh defined as 

_ 1 ^r-^ \P2num{ti) ~ P2model{ti)\ s 

with tiS are taken from t = 0.1 to t = 5 by 0.1 intervals. 

Figure [S] shows plots of rj versus interaction strength g plots for two different models. 
Model ( l48l) shows good agreements with simulation up to (? = 17, with relative error less 
than 3%. The error increases steadily, reaching more than 5% when g > 60 (seen from 
figure [21 this is strongly repulsive region). Second model f H9|) shows large deviation from 
the numerical simulation for smaller g, but agreements with the simulation becomes better 
than that of ( l48l) model for g > 67. 

Finally, The probability of both particles outside, P3, is easily calculated since Pi + P2 + 
P3 = 1. So total decay mechanism can be described by ( l43l) . ( HHl) or ( H9l) . 




1234567 ' '''''i''''2''''3''''4''''5''''6''''7 ' '''''i''''2''''3''''4''''5''''6''''7 * 



(a) g=0 (b) g=l (c) g=10 

FIG. 5: P2{t) plots from pHj) (line) and numerical simulation (dots) for g = 0,1, 10. Dots 
represent P2{t) from numerical simulation, dashed red lines are from f l48|) with 723 is 
substituted by lowest decay rate of single particle and blue lines are from (HHj) with 723 
obtained from fitting. All three show very good agreements and lines are almost 

overlapping. 
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(a) 



(b) 



FIG. 6: (a) The fitted decay rate 723 (solid line) and the lowest single particle decay rate 
7zO (dashed line) vs g. (b) the relative difference (723 — '^zo)hzo- 




1 2 3 4 5 6 

(b) g=100 



FIG. 7: Comparison between numerical simulation (dots) and theoretical models. Black 
line represents fl49p model and blue line represents (HSj) model. (H^ model shows better 
agreement with numerical simulation for larger g. 




20 40 60 80 100 120 140 

FIG. 8: Absolute mean of relative error 77 plots of model ( HHl) (blue) and model ( H9l) (thick 

black) versus g. 
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FIG. 9: InA^in(t) plots. For smaller g it is closer to the straight line (exponential decay) 
but for larger g the decay rate changes from faster to slower ones. 

If we calculate the number density Ni^{t) inside the potential, it is written as 

Nin{t) = dx dXidX2ij* {Xi,X2, t)'^^S{x — Xi)lp{Xi,X2,t) (52) 
Jin J 

and in our case it simply becomes Ninit) = 2Pi(t) + P2(t). Figure [9] Shows the logarithm 
of A^in versus time. For g = 0,1 the decay of Nin(t) is close to exponential (In A^in close to 
straight line) but for larger g it is more visible that the decay rate changes from faster to 
slower ones, as expected. 



V. CONCLUSION 



We have calculated the decay of two repulsively interacting bosons, initially in the ground 
state of a potential trap, by numerical simulation. We have found an exponential decay 
mode for the probability of both bosons inside the trap and have estimated its decay rate 
theoretically. By applying outgoing boundary condition for the loss of single particle from 
the trap, we obtain two complex wavevectors corresponding to the two particles inside 
the trap and corresponding decay rate. The agreement between numerical simulation and 
theoretical estimation in time evolution of decay probabilities is very good. For describing 
the probability to have a single particle inside and one particle outside, two simple models 
were proposed. For small and moderate g, we apply a model in which the remaining particle 
decays exponentially, whereas for larger g (strongly repulsive) we propose another model in 
which the modes of each complex wave vector decay separately. Our numerical simulations 
show very good agreement for weak and moderate interactions with the first model. For 
stronger interactions, where fermionization effects become relevant, separate exponential 
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decay model becomes necessary and agrees well with simulations. The number density shows 
that the decay rate changes over time from fast to slower decay for large g. The results show 
that simple models describe the overall decay mechanism of repulsively interacting bosons 
well. 
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